Self-potential time series reveal emergent behavior in soil organic matter dynamics

The active cycling of carbon between soil organic matter and the atmosphere is of critical importance to global climate change. An extensive body of research exists documenting the capricious nature of soil organic matter (SOM) dynamics, which is symptomatic of an intricate network of interactions between diverse groups of heterotrophic microorganisms, complex organic substrates, and highly variable local environmental conditions. These attributes are consistent with elements of complex system theory and the temporal evolution of otherwise unpredictable patterns of behavior that emerge from long range dependency on initial conditions. Here we show that vertical depth profile of self-potential (SP) time series measurements responds in a quantitative manner to variations in soil moisture, SOM concentrations, and relative rates of microbial activity. Application of detrended fluctuation analysis (DFA) of self potential time series data is shown additionally to reveal the presence of long-range dependence and emergence of anomalous electrochemical diffusion behavior, both of which diminish with depth as SOM specific energy densities decline.

The vast quantity of carbon stored as organic matter in the top 3 m of soils is greater than amounts that exist in the atmosphere and terrestrial vegetation combined [1][2][3][4] . Because of this, even a small change in the biogeochemical processes that contribute to the conversion of soil organic matter (SOM) into greenhouse gases such as carbon dioxide or methane could further exacerbate global climate change even as mean temperatures rise and extreme variations in rainfall become more frequent 5 . Higher turnover rates of SOM in surface soils and emission of greenhouse gases also threaten the vertical transport and long-term sequestration potential of organic carbon in deeper soil horizons where radiocarbon ages are commonly as old as one thousand to more than ten thousand years 1,4,6 . Investigation of these interconnected vulnerabilities in SOM dynamics has exposed significant challenges to researchers not only in terms of how to integrate traditional physical, chemical, and biological measurements in a quantitative manner, but also how to reduce the tremendous amount of time spent on sampling for multiple analyses. An alternative approach in soil hydrological and biogeochemical studies is the growing use of self-potential (SP) monitoring, a geophysical method based on passive measurement of electrical potential differences generated by the movement of charged ions and water through porous subsurface materials 7,8 .
SP has been used to locate ore deposits, map hydrothermal zones, monitor groundwater flow and characterize contaminated sites due to their associated electrochemical, thermoelectric, electrokinetic, piezoelectric, and redox effects 9 . The method holds considerable promise to be applied in new ways as demonstrated in this study, especially with respect to a time-series analysis of its signals which opens new opportunities for quantifying moisture dependent organic matter-microbial activity dynamics. In soils, natural electric potential gradients can be attributed to electrokinetic effects caused by flowing water and biochemical effects attributable to oxido-reduction phenomena 10 . Decoupling the source of SP signals is non-trivial but petrophysical relationship between SP signals and hydraulic gradient is well established through an electrokinetic coupling coefficient 9,11 . This approach facilitates signal interpretation by taking the coupled effects of mobile charged solutes, as well as the charged nature of stationary minerals and associated organic matter, into account 12 . SP signals therefore offer an indirect but robust approach to quantify soil processes including soil moisture, soil organic matter content and relative microbial activity driving biogeochemical changes in soils. In this study, we use strong correlations between measured SP signals and soil moisture, SOM, and relative microbial activity rates to demonstrate the capability of using SP to scale understanding of SOM dynamics in a non-invasive way with minimal disturbances. We also use a time-series analysis of SP data to demonstrate the transient behaviour of SOM-microbial activity rates.

Results and discussion
Soil hydraulic conditions and microbial activity. The moderately high saturated hydraulic conductivity values determined for Meilleurs Bay terrace sands are conducive to rapid infiltration and downwards movement of meteoric precipitation through the unsaturated zone (Fig. 1A). While depth variations in hydraulic conductivities parallel corresponding fining and coarsening trends in grain size (Fig. 1A), the well-sorted nature of the sands constrain porosity to a narrow range near 39% (Fig. 1B). Measurements of bulk density follow the same trends, decreasing in the direction of fining upwards grain sizes and increasing with coarsening upwards grain sizes, with the lowest bulk densities close to the surface (1482 kg/m 3 ) and 1.45 m underground (1443 kg/m 3 ). Considering the nearly constant porosity and uniform mineral composition of the terrace sands, the observed variations in bulk density can be attributed directly to differences in soil organic matter content with the greatest amounts predicted for depths corresponding to the lowest bulk densities (Fig. 1C) [13][14][15][16] . The peak of SOM at 1.45 m depth is consistent with the presence of a root mat horizon in the upwards coarsening interval between 0.9 and 1.5 m depth 17 .
Post-wetting drainage dynamics of unsaturated porous media are complex owing to spatial and temporal changes in the extent of hydraulic continuity between pores 12,18,19 . At the time of this investigation, 24 h after a rainfall event, moisture levels were close to the residual volumetric water content anticipated for rapid drainage of meteoric precipitation through highly permeable sandy soils (Fig. 1B) 18 . Elevated water contents at 0.10 and 1.45 m depths correspond directly to depth levels with the smallest grain size, lowest hydraulic conductivity, and highest content of soil organic matter. Among these parameters, linear regression established that soil organic matter content accounted for the greatest proportion of the variance in volumetric water contents (R 2 = 0.81; Fig. 1D). Similar relationships have been reported in other studies documenting the enhancement of water retention by soil organic matter 16,20 , including extracellular polymeric substances produced in situ by microorganisms 21 .
Decomposition of organic matter in soils by heterotrophic microbial respiration is strongly dependent on soil moisture 22,23 . This relationship is reflected in the similar depth profiles of calculated soil moisture microbial activity 24 and measured volumetric water contents ( Fig. 1E and B, respectively). Another underlying premise of the microbial activity-soil moisture function is that relative soil moisture microbial activity is reflective of the bioavailability of soil organic matter under partially saturated conditions 24 . This invokes the Monod kinetics model of microbial growth and substrate utilization under nutrient limiting conditions 22,25 . Recognizing that soil moisture and organic matter content are related (Fig. 1D) provides a useful paradigmatic link between the microbial activity-soil moisture function and bioavailability of organic matter implied by the Monod relationship. As shown in Fig. 1F, soil moisture microbial activity at different horizons across the depth profile is, in fact, responsive to organic matter concentrations, calculated independently as a function of bulk density. The implication is that soil moisture does limit the bioavailability of soil organic matter under fluctuating soil moisture conditions, with a half-saturation constant K = 0.97 + 0.03 for a relative soil moisture microbial activity equivalent to one-half of the level of microbial activity expected under optimal soil moisture conditions.  Fig. 2A). The trend reversal in SP response with depth is a mirror reflection of that observed with soil water content, soil organic matter, and soil moisture microbial activity. Of these parameters, water content accounts for 83% of the depth profile variation in SP (Fig. 2B) compared to 96% for soil organic matter (Fig. 2C) and 90% for soil moisture microbial activity (Fig. 2D).
The linear correlation between mean SP response and water content as an independent variable is consistent with numerous theoretical and experimental studies 9,11,12 ; however, considerable variation exists in the reporting of coupling coefficients (i.e., slope of the regression), depending on the saturated streaming potential coupling coefficient ( C s ) and how differences in pore water specific energy are expressed 10,26,27 . Here we apply a version of the Helmholtz-Smoluchowski equation for C s that incorporates the effective conductivity of a saturated porous medium, as defined by Archie's Law 7,12,28,29 ; relevant equations and parameters and are defined in the methods section. For an unsaturated porous medium, the difference in specific energy can be formulated in terms of volumetric water content as dP = L c ρ w gφ −1 dθ where L c is a soil-specific characteristic length scale that defines the maximum extent of hydraulic continuity for gravity drainage 30 ; C s L c ρ w gφ −1 thus represents a conditional unsaturated streaming potential coefficient. For the field data, this approach gave -0.024 mV/Pa for C s , a soil zeta potential of − 46.22 mV, and 0.10 for the volumetric water content ( θ ref ) at the shallow depth of the reference electrode. The first two values are consistent with those reported for fine quartz sands at low electrical conductivities (ca. 10 -3 S/m) that accompany infiltration of dilute meteoric precipitation 26,31 , while the estimate for θ ref is comparable to that measured at a depth of 0.10 m in the borehole (Fig. 1B).
The strong effect of SOM on SP is an intriguing finding that is related, in electrochemical terms, to surface charge development and cation exchange capacity of organic matter-mineral aggregates in soils, including those containing microbial biomass 12,[32][33][34] . Recognition of the relationship between SOM and volumetric water content (Fig. 1D) allows for refinement of the streaming potential coupling coefficient and extends the theoretical basis for quantitative interpretation of variations in SP response to observed differences in soil organic matter concentrations (Fig. 2C). In this case, the difference in specific energy takes the form dP = L c ρ w gm SOM φ −1 dSOM with m SOM equivalent to the slope of the linear regression of θ as a function of SOM; L c is a soil-specific length scale that defines the maximum extent of hydraulic continuity in the porous network as a function of parameters α and n of the van Genuchten model of soil water retention 30 . The corresponding estimates for C s and soil zeta potential from the SP data are − 0.028 mV/Pa and − 55.2 mV, respectively. These estimates compare favorably to those derived on the basis of SP response as a function of volumetric water content. At the same time, a calculated reference  www.nature.com/scientificreports/ electrode SOM concentration of 1.6% w/w is comparable to that at 0.10 m (1.53% w/w) and 1.45 m (1.41% w/w) depths in the borehole (Fig. 1C). As expected for this condition , the observed SP response at these depths is near 0 mV ( Fig. 2A). This illustrates the potential to apply inverse modeling in SP surveys to map out SOM concentrations in soils by geophysical means. The correlation between SP and soil moisture microbial activity presumably reflects the impact of changes in water content on mass transport processes and the manifest bioavailability of soil organic matter. Because optimal microbial activity is realized at a water saturation of 0.65 24,35 , the soil-specific characteristic length scale that defines the upper drainage limit at f SMA = 1.0 is taken to be 0.65L c 30 . This gives the difference in pore water specific energy in relation to microbial activity as dP = 0.65L c ρ w gdf h , which yields estimates of -0.027 mV and -53.2 mV for C s and zeta potential, respectively. A predicted value of f h = 0.62 was ascertained for soil moisture microbial activity at the reference electrode, in good agreement those calculated from measured water contents at 0.10 m ( f h = 0.60) and 1.45 m ( f h = 0.66) depths in the borehole (Fig. 1E). Accounting for SP response as a function of soil moisture microbial activity through an unsaturated streaming potential coupling coefficient provides a new quantitative geophysical perspective on how changes in pore water energy status are apt to influence the bioavailability and fate of SOM in soils.
Self-potential and long-range correlation. The application of DFA to examine the correlation structure of time series data has become a widely used computational method in time series analysis 17,36-38 . The SP time series yielded DFA α scaling exponents that declined with depth from as high as 1.75 in the upper 1.45 m of the profile to 0.90 at 3.5 m depth (Fig. 3A). This trend reflects a progressive weakening of long-range correla- The most important nonequilibrium driver of energy flow in soil systems is organic matter. While SOM exhibits a high degree of heterogeneity with respect to molecular structure and biogeochemical reactivity, specific energy density (J/mg SOM) from differential scanning and bomb calorimetry measurements provide a direct estimate of the energy state of SOM with respect to susceptibility of oxidation to carbon dioxide. SOM specific energy density can be modeled as decreasing exponential function of depth based on measurements across a wide variety of soil types [42][43][44][45][46][47] (Fig. 3A). This permits comparison of the SP time series DFA values with corresponding estimates of SOM energy densities (Fig. 3B). The trend evident in the data follows a distinctive pattern of a logistic transition between energy states with α values increasing as a function of SOM energy density 48 . The inflection point of the regression curve at -0.59 J/mg SOM is at an energy level where the predicted range of DFA α values evoke the onset of Brownian motion and normal diffusion in the SP fluctuations. Conversely, high SOM energy densities and α scaling exponents fall into the range of super-diffusive electrochemical behavior that is symptomatic of contributions from high rates of microbial respiration 37,49 . Moreover, a cross-over from super-diffusion ( α > 1.0) to sub-diffusion ( α < 1.0) is apparent across a narrow range of SOM energy densities (− 2.0 to − 4.0 J/mg SOM) that implies a transition from non-equilibrium microbial catalyzed to mass action controlled chemical reactivity (e.g., mineral-SOM interactions) as the primary determinant of electrochemical energy flow in deeper subsoil horizons 48,49 .

Emergent behavior in soil organic matter dynamics.
There is ample evidence supportive to extending the paradigm of complex system theory to SOM dynamics 3,32,50 . Recognition of a high degree of spatiotemporal variation of interactions between physical, chemical, and microbial processes under non-equilibrium conditions in soils is an especially compelling consideration. Another hallmark of complex systems is the manifestation of emergent behavior in the form of a perceptible degree of long-range dependence and sensitivity to initial conditions within a rich network of interacting processes 5,6,51,52 . These attributes are mirrored not only by the observed correlations of self potential to soil moisture, SOM content and microbial activity (Fig. 2), but also the detection of long-range dependency and emergence of anomalous diffusion behavior in the SP time series data. Moreover, the progressive loss of long-range dependence with depth and approach to normal diffusion behavior is consistent with a weakening and possible loss of some network interactions with declining SOM energy densities. www.nature.com/scientificreports/ Another insightful aspect of the logistic relationship between DFA α values of the SP time series and SOM energy density is the changing width of the nonlinear regression prediction interval (Fig. 3B). Of particular interest is the greater range of the prediction interval ( α ∈ [1.45, 1.90]) at energy densities ≤ -3.0 J/mg SOM, which is representative of shallow depths ≤ 1.20 m. Considering that α depicts the strength of long-range dependence and emergence of anomalous diffusivity, the wide prediction interval suggests that potentially small perturbations in the underlying network of interacting dynamic variables (e.g., soil moisture, SOM content, microbial activity) could trigger almost an order of magnitude shift in emergent behavior. The potential for such unpredictable "chaotic" changes is perceived as another characteristic feature of temporal dynamics in complex systems, including heterogeneous soil environments 5,6,53 . Conversely, the much-narrowed prediction interval at SOM energy densities ≥ − 2.0 J/mg SOM signifies a weakening of long-range dependency and dissipation of anomalous diffusion behavior owing to a fragmentation of network interactions.
Conclusion. The correlation between SP and SOM, soil moisture and soil microbial activity demonstrated in this study highlights a new opportunity to quantify soil organic matter and soil moisture-microbial activity. The complexity with direct measurements of these soil properties makes alternative approaches valuable for constraining and scaling up the spatial and temporal resolution of their estimates. While this study uses borehole SP measurements, surface measurement of SP has been demonstrated for groundwater and contaminated sites characterization. Hence, SP provides a non to minimally invasive geophysical tool which when combined with minimal direct measurement would improve the characterization and monitoring the spatial and temporal variations in soil organic matter and microbial activity. This is bound to expand our mechanistic understanding of the capricious nature of soil organic matter. Besides elucidating the transient nature of electrochemical energy flow in SOM relatable to soil microbial respiration, the time-series SP data acquisition presented in the studies challenges the conventional approach of electrical geophysical data acquisition broadly applied in subsurface characterization and monitoring including complex environmental systems 9,10,12,38 . Specifically, electrical signatures are commonly measured multiple times and an average value calculated and attributed to a given location. In complex and emergent system, this ignores the pronounced temporality that may be inherent in such system. We therefore recommend adoption of a time-series measurement approach for soil electrical signals (mostly natural potentials) as this is apt to capture temporal perturbations indicative of complex time-varying emergent processes within the subsurface.

Methods
Study site description. The study site is a small forest clearing (NAD 83 UTM coordinates 18 T 297593E 5115217 N; elevation 133 masl) in the upper reaches of the Meilleurs Bay catchment on the south bank of the Ottawa River, approximately 10 km northwest of Deep River, Ontario, Canada 13 . The mean annual temperature in the region is 4.9° C, ranging from 18.4 °C in July to − 10.5 °C in January 54 . The average annual precipitation is 878 mm with the local climate falling into the Köppen classification scheme as Dfb, humid continental [55][56][57] . The forest at the study site is composed of balsam fir (Abies balsamea), red maple (Acer rubrum), yellow birch (Betula alleghaniensis), white spruce (Picea glauca), eastern white pine (Pinus strobas), and poplar (Populus tremuloides). Prominent species along the edge of the clearing are wild blueberry (Vaccinium angustifolium), black raspberry (Rubus occidentalis), and wild strawberry (Frangaria vesca). Ground vegetation in the clearing includes fescues and tufted grasses (Festuca and Lolium sp.), goldenrod (Solidago sp.), sweet fern (Comptonia peregrina) and hawkweed (Hieracium sp.).
Around 10.5 ka BP, the retreat of the Laurentide Icesheet from the Precambrian Canadian Shield opened the Ottawa River valley to drainage from pro-glacial meltwater lakes in the Huron and Michigan sub-basins of the Great Lakes 13,58 . This discharge contributed to extensive deposition of medium-to fine-grained glaciofluvial channel sands along the slopes of the river valley. As water levels fell in response to incremental icesheet recession and changes in drainage routes of the pro-glacial lakes, subaerial exposure of the channel sand deposits contributed to the development of humo-ferrric podzol soils (Canadian soil classification scheme), which dominate over much of the upper Ottawa Valley region, including the Meilleurs Bay locality 4,17 . At the Meilleurs Bay, a series of terraces between 129 and 141 masl document brief fluctuations in river height during the recession of the Laurentide Ice Sheet 36 .
Instrumentation of the Meilleurs Bay catchment with an extensive network of piezometers has established that the unsaturated vadose zone at the study site extends to a depth of 5.0 m or more 59 . Seasonal and interannual fluctuations in water table depths throughout the aquifer are small, typically less than ± 0.5 m. Regional hydrological studies indicate that approximately 37 percent of the annual precipitation contributes to groundwater recharge or runoff. The remainder is returned directly to the atmosphere by evapotranspiration. Almost all the groundwater recharge occurs during and just after the spring snowmelt, whereas strong evapotranspiration during the summer months gives rise to soil moisture deficits.

Soil depth profile and self-potential measurements.
An open-barrel stainless steel sand auger with a diameter of 5 cm was used to advance a borehole to a final depth of 3.5 m through the unsaturated vadose zone at the study site. Sediment samples were collected for determination of grain size distribution, bulk density, and volumetric water content at 0.15 to 0.50 m intervals. Sample temperatures were recorded immediately after collection using a Solinst Levelogger to facilitate estimation of prevailing water density ( ρ w ) and dynamic viscosity ( µ ) values. Saturated hydraulic conductivities were calculated from the grain size distributions using a modified Hazen equation 60 www.nature.com/scientificreports/ with the coefficient of uniformity ( C u ) equal to the ratio between the 60th ( d 60 ) and 10th ( d 10 ) percentile particle size diameters (in meters) and g corresponding to the acceleration of gravity. Porosity ( φ ) was calculated as a function of C u 24,47 .
Bulk densities ( ρ b ) were used to estimate soil organic matter concentrations (SOM % by weight) following the relationship derived for sandy forest soils of glacial fluvial origin 10,26 Measurements of self potential were made at each depth interval using a Fluke 289 logging multimeter with an internal impedance of 20 MΩ. A pair of Cu-CuSO 4 non-polarizing electrodes were connected to the terminals of the multimeter by two 30 m lengths of 18 gauge solid-core insulated copper wire. The lead electrode was fastened to a length of 2.5 cm outer diameter PVC conduit (to permit installation and retrieval from the borehole) and connected to the positive terminal of the multimeter. The reference electrode was connected to the negative terminal and inserted 0.15 m into the ground approximately 40 m away from the borehole. Readings were recorded as a time series at a frequency of 1 Hz over an interval of 1200 s (20 min). Potential differences recorded by the multimeter represent an average of five readings taken over the course of one second.
Assessment of microbial activity. Relative rates of microbial heterotrophic respiration were calculated as a function of volumetric water contents using the microscale soil moisture microbial activity function ( f h ) developed by Yan et al. 24 The value of f h compares heterotrophic respiration rates ( r θ ) at different levels of soil moisture ( θ ) to the maximum rate ( r max ) sustained by an optimal water content ( θ opt ), which is taken to be approximately 65 percent of soil porosity ( φ ) (i.e., θ opt = 0.65φ ). Suggested values of K θ = 0.1 for the soil organic matter desorption constant and n = 2 for the saturation exponent were used along with the soil organic matter-microorganism collocation factor a , estimated as a function of fine grain mineral content (< 45 µm diameter fraction; a = 2.8c 45 − 0.046 ), to complete parameterization of the f h function 24 .
At soil moisture values below 0.65φ , the loss of hydraulic continuity between pores restricts soil organic matter bioavailability and limits heterotrophic microbial activity as predicted by the Monod microbial growth model 24,25 .

Rearrangement yields
Here γ is the yield coefficient for microbial growth on soil organic matter, µ max is the maximum microbial growth rate constant, m is the microbial biomass concentration, [SOM] is the soil organic matter concentration (%), estimated from bulk density values, and K S is the half-saturation constant (i.e., the SOM concentration where the rate of heterotrophic respiration is 0.5 r max ).
Evaluation of self-potential response. In a saturated porous medium, the difference in self-potential ϕ with respect to specific energy of pore water, expressed in terms of pressure ( P , energy per unit volume) is proportional to a streaming potential coupling coefficient C s The Helmholtz-Smoluchowski equation gives C s as a function of the electrical permittivity of water ( ǫ ), zeta potential ( ζ ), dynamic viscosity ( µ ), and effective conductivity ( σ) The influence of pore space connectivity on effective conductivity is described by Archie's Law www.nature.com/scientificreports/ where σ w is the electrical conductivity of the water, φ is porosity, and m is the porosity (or cementation) exponent 12,28 . This gives As water saturation ( S w = θ/φ ) decreases during drainage, the change in pore water specific energy is characterized by a soil-specific length scale ( L c ) that defines the maximum extent of hydraulic continuity in the porous network as a function of parameters α and n of the van Genuchten model of soil water retention Within the limits of drainage between saturated and residual water contents of a soil Rearrangement and substitution into (7) yields which gives a linear relationship for self-potential differences as a function of volumetric water contents at reference ( θ ref ) and measurement ( θ ) electrodes.
Energy density of soil organic matter. The energy density ( E d ) of soil organic matter tends to decrease in an exponential manner with respect to depth z 43 , with a decay constant k = 1.35 + 0.19 m-1, and surface E d(z=0) = − 17.70 + 0.60 J/mg SOM ascertained for 125 measurements of SOM E d at different depths across a wide variety of soil types [42][43][44][45][46][47] . This empirical relationship was applied to derive estimates for SOM E d values at borehole depths corresponding to each of the measured SP time series. Data analyses. Nonlinear estimation using Levenberg-Marquardt optimization in STATISTICA 13.2 was employed to derive parameter estimates and the portion of variance accounted ( R 2 ) for in dependent variables of regression equations. Analyses of the self-potential time series were performed in Python (3.8) by detrended fluctuation analysis (DFA) using the dynamical systems package Nolds (0.5.2) following the procedure described by Peng et al. 27 ; the minimum time interval of windows was limited to 60 s in length, and computations were performed using both non-overlaping windows and successsive windows overlapped by half of their length. To further refine our estimate of DFA precision, each time series was downsampled into two time series equivalent to a sampling frequency of 0.5 Hz and the scaling exponents of the two new time series were calculated independently as described above. For all of the time series, the standard error of mean DFA α values was < ± 0.06. Random shuffling of each time series yielded α = 0.5, as expected for an uncorrelated signal.

Data availability
Data available on request from the corresponding author. www.nature.com/scientificreports/